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A new algorithm for numerical integration of the rigid-body equations of motion is proposed. The 
algorithm uses the leapfrog scheme and the quantities involved are angular velocities and orienta- 
tional variables which can be expressed in terms of either principal axes or quaternions. Due to 
specific features of the algorithm, orthonormality and unit norms of the orientational variables are 
integrals of motion, despite an approximate character of the produced trajectories. It is shown that 
the method presented appears to be the most efficient among all known algorithms of such a kind. 

PACS numbers: 02.60.Cb; 95.75.Pq; 04.25.-g 
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The method of molecular dynamics (MD) plays a 
prominent role in studying molecular liquids. All exist- 
ing techniques appropriate to simulate such systems can 
be categorized in dependence on what type of parameters 
are chosen to represent the rotational degrees of freedom 
and what kind of numerical algorithm is applied to inte- 
grate the corresponding equations of motion. 

In the molecular approach, the phase trajectories are 
considered in view of translational and rotational mo- 
tions. The translational dynamics is defined by motion 
of molecular centers of masses, whereas the orientational 
motion can be determined in terms of Eulerian angles [1, 
2], quaternions [3-8] or principal-axis vectors [4]. The 
numerical integration within Eulerian angles is very in- 
efficient due to singularities of the equations of motion 
[3, 5]. If the quaternions or principal- axis vectors are in- 
volved, additional efforts must be paid to conserve their 
unit norms or orthonormality. 

The atomic approach [9] treats dynamics of the sys- 
tem in view of translational motion of individual atoms 
which move under the potential-energy forces plus forces 
of constraints introduced to hold inter-atomic distances 
constant. This approach is believed to have good stabil- 
ity properties, because the usual Verlet algorithm can be 
applied here. Nevertheless, the atomic approach is so- 
phisticated to implement for point molecules and when 
there are more than two, three or four atoms in the cases 
of linear, planar and three-dimensional molecules, respec- 
tively. Moreover, to reproduce the rigid molecular struc- 
ture it is necessary to solve complicated systems of non- 
linear (in general, six per molecule) equations at each 
time step of the integration [10]. 

It is a common practice to integrate orientational mo- 
tion with the Gear predictor-corrector algorithm of a 
high-order [11]. Such an algorithm, being accurate at 
very small time steps, quickly becomes unstable with in- 
creasing the step size [10]. Translational motion is usu- 
ally integrated with lower-order Verlet [12], velocity Ver- 
let [13] and leapfrog [14] integrators, owing their simplic- 
ity and exceptional numerical stability. However, original 
versions of these integrators were constructed assuming 



that acceleration is velocity-independent and, therefore, 
they can not be applied directly to rotational dynamics. 
Analogous problems arise with translational motion in 
the presence of magnetic fields. 

In order to remedy that situation, Fincham [15] has 
derived a rotational-motion version of the leapfrog algo- 
rithm in which systems of four nonlinear equations per 
molecule for quaternion components are solved by itera- 
tion. Ahlrichs and Brode have introduced a method [16] 
in which principal axes are considered as pseudo-particles 
and constraint forces are introduced to maintain their 
orthonormality. But the algorithm is within the Verlet 
framework and does not contain angular velocities ex- 
plicitly. The quaternion dynamics with constraints was 
also formulated [17]. As a result, a new algorithm within 
the velocity Verlet framework has been generated. Re- 
cently, the principal-axes scheme has been adapted to 
this framework as well [18]. Nevertheless, it was con- 
cluded that the best numerical stability can be achieved 
in the atomic-constraint approach. 

In this paper we propose a new leapfrog integrator of 
the rigid-body equations of motion. The main idea con- 
sists in involving angular velocities, instead of angular 
momenta, into the integration. This leads to significant 
simplifications with respect to angular-momenta versions 
[15]. The algorithm seems to be the most efficient and 
simple, exhibiting excellent stability properties which are 
similar to those observed within the cumbersome atomic- 
constraint technique. 

Consider a classical system with N rigid molecules 
composed of M point atoms. Translational motion of 
the system is described in the usual way, applying New- 
ton's law, whereas two first-order equations per molecule 
of the rotational dynamics can be obtained as follows. 
According to Eulcr equations [1], the rate of change in 
time of principal components, (fix> ^V> ^z) = °f an ~ 
gular velocity is 
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= K i a (t) + (jp-j^n^(t)sii y (t). 
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Here (a,/3, 7 ) = (X,Y,Z), (Y,Z,X) and (Z,X,Y), K\ 



1 



are principal components, Kj = A^k^, of the torque 

kj = Y^^-'^l ( r i ~ r i) x ^ij' exerted on molecule i with re- 
spect to its center of mass due to the site-site interac- 
tions ffj = f (r° — Tj) with the other molecules, J a denote 
the principal moments of inertia, oricntational variables 
were collected into the square orthonormal matrices Aj, 
the nine elements of each of which (i = 1, . . . , N) present 
coordinates of three principal axes (XYZ) of the molecule 
in the laboratory frame, the position of atom a within 
molecule i in the same frame is rf(t) = r;(t) + A+ (i)A a , 
where A a = (A", A£, Ag) + is a vector-column of these 
positions in the body frame and A + the matrix trans- 
posed to A. 

The second equation follows from definition of angular 
velocity, 



dAi 
dt 
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Aj=W(ftj)Aj, (2) 



where the property AA + = I of orthonormal matrices 
has been used, W(fij) is a skewsymmetric matrix asso- 
ciated with angular velocity, i.e., W + (f2j) = — W(flj) 
and I designates the unit matrix. In an alternative rep- 
resentation the matrix Aj = A(q^) is a function of the 
four-component quaternion qi = (£j, rji, Q, Xi) + [4, 5]. 
The time derivatives of quaternions can be cast in the 
form 
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q t = Q(fij)q i; 



(3) 



where Q(fij) is a skewsymmetric matrix again and the 
unit quaternion norm + rjf + Q + xf = 1 , which follows 
from the orthonormality of Aj, has been used. 

In the case of translational motion, it is easy to derive 
the leapfrog algorithm [14]: Vj(t+|) = Vj(i- f ) + /iaj(t), 
Yi(t + h) = Yi(t) + hviit + -|), where h denotes the time 
increment, v, = dri/dt is the center-of-mass velocity, 
a.i(t) = ■^Y^j-'c^l^ijit) tnc molecular acceleration and 
m the mass of a separate molecule. Recently, it has been 
shown that contrary to the conventional point of view, 
the order of truncation errors for this leapfrog is four 
rather than three for both coordinates and velocities due 
to a fortunate cancellation of uncertainties [19]. 

The problems with deriving a leapfrog algorithm for 
rotational motion are that angular accelerations (1) de- 
pend explicitly not only on spatial coordinates via molec- 
ular toques but also on angular velocities. Moreover, the 
time derivatives of orientational variables do not define 
angular velocities directly (see Eqs. (2) and (3)). These 
difficulties can not be handled with a simple leapfrog in 
which position and velocity are known at different times. 
It is worth to underline that similar problems (even much 



more difficult) arise in the angular-momentum approach 
[15], Verlet and velocity Verlet frameworks [17, 18]. 

The basic idea of our approach lies in involving princi- 
pal angular velocities into the integration process. Then, 
acting in the spirit of leapfrog scheme and using Euler 
equation (1), one obtains 



+ (J P - J 7 )fij, (n) (t)fif°(t) 



(4) 



While the molecular torques K l a {t) can easily be evalu- 
ated via the coordinates rj {€) and Aj (t) or qi (t) , a prop- 
agation of the products of angular velocities in Eq. (4) 
to on-step levels of time is necessary. The obvious choice 
for this is 



ni (n) (t)fi* w (t) = ^ ni(t-i)ni(t-§) 
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(5) 



In view of (5), equation (4) constitutes a system of max- 
imum three nonlinear equations per molecule for the un- 
knowns fl l a (t+ 1 ) . The system is simple and can be solved 
in a quite efficient way by iteration, n = 0, 1, . . . , taking 

^a \t + f ) = fi a(* - |) as an initial guess. We note 
that the order of truncation errors for angular-velocity 
evaluation (4) reduces to three, because approximation 
(5) is only second order accurate on h. 

The evaluation of orientational variables can be real- 
ized by writing 



Si{t + h) = S i (t) + hH i S i {t+§) 



(6) 



for principal-axis vectors (Si = Aj,Hj = Wj) and 
quaternion (Si = qi,Hi = Qi) representations, where 
Eqs. (2) and (3) have been used. The matrices Wj = 
W(fij) and Qi = Q(fi,) are calculated using already de- 
fined angular velocities Jlj(i + -|), whereas orientational 
variables can be propagated to mid-step levels of time as 



S*(* + #) = ±[Si{t) + Si(t + h)]. 



(7) 



Equation (6) together with (7) are, in fact, systems of lin- 
ear equations with respect to elements of Aj(t + h) and 
qi(t + h), which, therefore, can be solved analytically. 
The result is 

Sj(t + h) = (I - f Hj)" 1 ^ + |Hi)Si(i) = @i{t, h)Si(t). 

(8) 

More explicit expressions for the set 0i = {Di,Gi} of 
evolution matrices are: Di = [1(1 — ^-fif) + hWi + 

£ Pj ]/[l+i£n?] andGi = [I(l-f^?)+/ ! Q,]/[l + ^ 2 ] 
in the cases of principal axes and quaternion representa- 
tions, respectively, where Pi is a symmetric matrix with 
the elements WJt + §)fij,(i +§) and fi? = «?(*+ §). 



2 



This completes the algorithm. It is interesting to remark 
that evaluation (8) exhibits the same fourth-order local 
accuracy on h as in the case of translational coordinates, 
despite the second order of interpolation (7). The reason 
for this results again from a cancellation of errors arising 
in coordinates and velocities during two neighbor time 
steps. 

It can be verified easily that the matrix (I — AH) _1 (I + 
AH) is orthonormal at arbitrary values of A, provided 
H+ = H. Then, as follows from construction (8), the 
evolution matrices and Gj are orthonormal as well. 
Therefore, if initially the orthonormality of A, and unit 
norms of are satisfied, they will be fulfilled perfectly at 
arbitrary times in future, despite the approximate char- 
acter of produced trajectories. This fact can be consid- 
ered as the main advantage of the algorithm derived that 
distinguishes it from all other singularity free algorithms, 
because no additional efforts are needed to preserve the 
rigid structure of molecules. 

We now test our approach on the basis of MD sim- 
ulations on liquid water. The simulations were per- 
formed in an NVE ensemble with N = 256 molecules 
at a density of N/V—l g/cm 3 and at a temperature of 
298 K using the TIP4P potential (M = 4) and reac- 
tion field geometry [20]. All runs were started from an 
identical well equilibrated configuration. The numeri- 
cal stability was identified in terms of fluctuations of 
the total energy, £ = \{{E - {E)) 2 )] 1 ' 2 /\{E)\. The ki- 
netic part of the energy was calculated at time t putting 
V (t) = l[V(t-l)+V(t+l)} + 0(h 2 ) for V = {v^fi,}, 
where the main term 0(h 2 ) of uncertainties is in the self- 
consistency with the second order of global errors for our 
algorithm (one order lower than minimal order of trun- 
cation errors for coordinates and velocities). 

As the atomic-constraint algorithm [9, 10] is inten- 
sively exploited and its performances are generally recog- 
nized, we have made comparative tests using this method 
and our advanced leapfrog algorithm within quaternion 
and principal-axes variables, as well as all known other 
approaches, namely, the fifth-order Gear algorithm [11], 
implicit leapfrog of Fincham [15], pseudo-particle for- 
malism [16], quaternion- and matrix-constraint methods 
[17,18]. Samples of £(t) as a function of the length of 
the simulations at four fixed values of ft = 1, 2, 3 and 4 
fs are shown in Fig. 1 . The usual value of step size for 
studying such a system is 2 fs [21]. 

Despite the Gear algorithm integrates the equations of 
motion very well at h = 1 fs, it has a very small region 
of stability and can not be used for greater time steps 
(see Fig. 1 (b)). Small step sizes are impractical in cal- 
culations because too much expensive computer time is 
required to cover the sufficient phase space. At the same 
time, the quaternion- and matrix-constraint methods as 
well as the pseudo-particle approach produce much more 
stable trajectories and exhibit similar equivalence in the 
energy conservation. Worse results are observed for the 
Fincham's leapfrog method. Finally, the best numeri- 



cal stability is achieved in the atomic-constraint tech- 
nique and our leapfrog scheme within both quaternion 
and principal axes representations, which conserve the 
energy approximately with the same accuracy (the re- 
sults for principal-axis variables and pseudo-particle for- 
malism are not included in the figure to simplify graph 
presentation). Quite a few iterations (the mean num- 
ber of iterations varied from 3 to 5 at ft = 1 -j- 4 fs) 
was sufficient to find solutions to the system of nonlinear 
equations (4) with a precision of 10~ 12 . This contributes 
a negligible small computation time additionally into the 
total time. 

No shift of the total energy was observed for the 
atomic-constraint and our leapfrog techniques at ft < 4 
fs over a length of 10 000 steps. To reproduce features of 
an NVE ensemble quantitatively, it is necessary for the 
ratio r = £/T of the total energy fluctuations to the fluc- 
tuations T of the potential energy to be no more than a 
few per cent. We have obtained the following levels of £ 
at the end of the runs in our leapfrog approach: 0.0016, 
0.0065, 0.015 and 0.029 %, corresponding to T « 0.29, 
1.2, 2.7 and 5.2 % at ft,= 1, 2, 3 and 4 fs, respectively (for 
the system under consideration T « 0.56%). Therefore, 
the greatest time step considered (4 fs) is still suitable 
for precise calculations. The ratio T can be fitted with a 
great accuracy to the function Ch 2 with a coefficient of 
C ~ 0.29 % fs~ 2 . This is completely in line with our the- 
oretical prediction about a characteristic square growth 
of global errors and, as a consequence, £{t) at t 3> ft. The 
square growth was observed in all other approaches, ex- 
cepting the Gear algorithm. However, only the advanced 
leapfrog algorithm provides a minimum of C and total 
energy fluctuations. 

The algorithm presented might become popular be- 
cause of its great stability, simplicity to implement for 
arbitrary rigid bodies and its intrinsic conservation of 
rigid structures. These features should be considered as 
significant benefits of the algorithm with respect to all 
the rest approaches. It can easily be substituted into 
existing MD programs on rigid polyatomics. Moreover, 
since velocities appear explicitly, the algorithm can be 
extended to a thermostat version and to integration in 
the presence of magnetic fields. These problems will be 
discussed in a separate publication. 

The author thanks the President of Ukraine for finan- 
cial support. 
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Figure caption 

Fig. 1. The total energy fluctuations as functions of 
the length of the simulations on liquid water, performed 
in various techniques at four fixed time steps: (a) 1 fs, 
(b) 2 fs, (c) 3 fs and (d) 4 fs. 
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